######################################################################
## Business as Usual? Conventional Politics, Critical Junctures     ##
## and Policy Feedback in the Paycheck Protection Program           ##
##                                                                  ##
## Michael Schwan, Marc Schneiber and Mark Cassell                  ##
## May 29, 2025                                                     ##
##                                                                  ##
## NOTE                                                             ##
## This file contains the code for replicating our key regression   ##
## results as well as the accompanying coefficient plots:           ##
## Tables E1-E2 in the online appendix; figures 5-6 in the paper.   ##                                              ##
## Further details, code and data on the replication of additional  ##
## robustness checkes and graphs are available upon request.        ##
######################################################################

# Our main data set on the distribution of PPP loans across 
# US Congressional Districts is: POP_D_24_00450_REPLICATION

# We use the following R packages for analysis and
# visualization: plm, tidyverse, broom, broom.mixed,
# sandwich, stargazer, sjPlot, grid, gridExtra.

# First, we transform our main data set into a panel data frame
panel <- pdata.frame(pop_d_24_00450_replication,
                     index = c("CD", "period"))

# We use this panel to run the set of models in table E1
m1 <- plm(loan_rate ~ 
  period + log_covid + estab_covid_pct + log_pop + 
  rural_sub + sparse_sub + dense_sub + urban_sub + purban +
  state_AL + state_AR + state_AZ + state_CA + state_CO + 
  state_CT + state_DE + state_FL + state_GA + state_HI + 
  state_IA + state_ID + state_IL + state_IN + state_KS + 
  state_KY + state_LA + state_MA + state_MD + state_ME + 
  state_MI + state_MN + state_MO + state_MS + state_MT + 
  state_NC + state_ND + state_NE + state_NH + state_NJ + 
  state_NM + state_NV + state_NY + state_OH + state_OK + 
  state_OR + state_PA + state_RI + state_SC + state_SD + 
  state_TN + state_TX + state_UT + state_VA + state_VT +
  state_WA + state_WI + state_WV + state_WY,
  data = panel, model = "random")

m2 <- plm(loan_rate ~ 
  period + log_covid + estab_covid_pct + log_pop + 
  party + hwm + hsb + happ + ssb + pvi_abs +
  rural_sub + sparse_sub + dense_sub + urban_sub + purban +
  state_AL + state_AR + state_AZ + state_CA + state_CO + 
  state_CT + state_DE + state_FL + state_GA + state_HI + 
  state_IA + state_ID + state_IL + state_IN + state_KS + 
  state_KY + state_LA + state_MA + state_MD + state_ME + 
  state_MI + state_MN + state_MO + state_MS + state_MT + 
  state_NC + state_ND + state_NE + state_NH + state_NJ + 
  state_NM + state_NV + state_NY + state_OH + state_OK + 
  state_OR + state_PA + state_RI + state_SC + state_SD + 
  state_TN + state_TX + state_UT + state_VA + state_VT +
  state_WA + state_WI + state_WV + state_WY,
  data = panel, model = "random")

m3 <- plm(loan_rate ~ 
  period + log_covid + estab_covid_pct + log_pop + 
  log_don + log_estab_large +
  rural_sub + sparse_sub + dense_sub + urban_sub + purban +
  state_AL + state_AR + state_AZ + state_CA + state_CO + 
  state_CT + state_DE + state_FL + state_GA + state_HI + 
  state_IA + state_ID + state_IL + state_IN + state_KS + 
  state_KY + state_LA + state_MA + state_MD + state_ME + 
  state_MI + state_MN + state_MO + state_MS + state_MT + 
  state_NC + state_ND + state_NE + state_NH + state_NJ + 
  state_NM + state_NV + state_NY + state_OH + state_OK + 
  state_OR + state_PA + state_RI + state_SC + state_SD + 
  state_TN + state_TX + state_UT + state_VA + state_VT +
  state_WA + state_WI + state_WV + state_WY,
  data = panel, model = "random")

m4 <- plm(loan_rate ~ 
  period + log_covid + estab_covid_pct + log_pop + 
  black + hispanic + pov +
  rural_sub + sparse_sub + dense_sub + urban_sub + purban +
  state_AL + state_AR + state_AZ + state_CA + state_CO + 
  state_CT + state_DE + state_FL + state_GA + state_HI + 
  state_IA + state_ID + state_IL + state_IN + state_KS + 
  state_KY + state_LA + state_MA + state_MD + state_ME + 
  state_MI + state_MN + state_MO + state_MS + state_MT + 
  state_NC + state_ND + state_NE + state_NH + state_NJ + 
  state_NM + state_NV + state_NY + state_OH + state_OK + 
  state_OR + state_PA + state_RI + state_SC + state_SD + 
  state_TN + state_TX + state_UT + state_VA + state_VT +
  state_WA + state_WI + state_WV + state_WY,
  data = panel, model = "random")

m5 <- plm(loan_rate ~ 
  period + log_covid + estab_covid_pct + log_pop + 
  party + hwm + hsb + happ + ssb + pvi_abs +
  log_don + log_estab_large +
  rural_sub + sparse_sub + dense_sub + urban_sub + purban +
  state_AL + state_AR + state_AZ + state_CA + state_CO + 
  state_CT + state_DE + state_FL + state_GA + state_HI + 
  state_IA + state_ID + state_IL + state_IN + state_KS + 
  state_KY + state_LA + state_MA + state_MD + state_ME + 
  state_MI + state_MN + state_MO + state_MS + state_MT + 
  state_NC + state_ND + state_NE + state_NH + state_NJ + 
  state_NM + state_NV + state_NY + state_OH + state_OK + 
  state_OR + state_PA + state_RI + state_SC + state_SD + 
  state_TN + state_TX + state_UT + state_VA + state_VT +
  state_WA + state_WI + state_WV + state_WY,
  data = panel, model = "random")

m6 <- plm(loan_rate ~ 
  period + log_covid + estab_covid_pct + log_pop + 
  black + hispanic + pov +
  party + hwm + hsb + happ + ssb + pvi_abs +
  rural_sub + sparse_sub + dense_sub + urban_sub + purban +
  state_AL + state_AR + state_AZ + state_CA + state_CO + 
  state_CT + state_DE + state_FL + state_GA + state_HI + 
  state_IA + state_ID + state_IL + state_IN + state_KS + 
  state_KY + state_LA + state_MA + state_MD + state_ME + 
  state_MI + state_MN + state_MO + state_MS + state_MT + 
  state_NC + state_ND + state_NE + state_NH + state_NJ + 
  state_NM + state_NV + state_NY + state_OH + state_OK + 
  state_OR + state_PA + state_RI + state_SC + state_SD + 
  state_TN + state_TX + state_UT + state_VA + state_VT +
  state_WA + state_WI + state_WV + state_WY,
  data = panel, model = "random")

m7 <- plm(loan_rate ~ 
  period + log_covid + estab_covid_pct + log_pop + 
  black + hispanic + pov +
  log_don + log_estab_large +
  rural_sub + sparse_sub + dense_sub + urban_sub + purban +
  state_AL + state_AR + state_AZ + state_CA + state_CO + 
  state_CT + state_DE + state_FL + state_GA + state_HI + 
  state_IA + state_ID + state_IL + state_IN + state_KS + 
  state_KY + state_LA + state_MA + state_MD + state_ME + 
  state_MI + state_MN + state_MO + state_MS + state_MT + 
  state_NC + state_ND + state_NE + state_NH + state_NJ + 
  state_NM + state_NV + state_NY + state_OH + state_OK + 
  state_OR + state_PA + state_RI + state_SC + state_SD + 
  state_TN + state_TX + state_UT + state_VA + state_VT +
  state_WA + state_WI + state_WV + state_WY,
  data = panel, model = "random")

m8 <- plm(loan_rate ~ 
 period + log_covid + estab_covid_pct + log_pop + 
 black + hispanic + pov +
 party + hwm + hsb + happ + ssb + pvi_abs +
 log_don + log_estab_large +
 rural_sub + sparse_sub + dense_sub + urban_sub + purban +
 state_AL + state_AR + state_AZ + state_CA + state_CO + 
 state_CT + state_DE + state_FL + state_GA + state_HI + 
 state_IA + state_ID + state_IL + state_IN + state_KS + 
 state_KY + state_LA + state_MA + state_MD + state_ME + 
 state_MI + state_MN + state_MO + state_MS + state_MT + 
 state_NC + state_ND + state_NE + state_NH + state_NJ + 
 state_NM + state_NV + state_NY + state_OH + state_OK + 
 state_OR + state_PA + state_RI + state_SC + state_SD + 
 state_TN + state_TX + state_UT + state_VA + state_VT +
 state_WA + state_WI + state_WV + state_WY,
 data = panel, model = "random")


cov.m1.HC <- vcovHC(m1, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m2.HC <- vcovHC(m2, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m3.HC <- vcovHC(m3, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m4.HC <- vcovHC(m4, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m5.HC <- vcovHC(m5, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m6.HC <- vcovHC(m6, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m7.HC <- vcovHC(m7, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m8.HC <- vcovHC(m8, method = "white1", type = "HC3", 
                    cluster = "group")

rob.se.m1.HC <- sqrt(diag(cov.m1.HC))
rob.se.m2.HC <- sqrt(diag(cov.m2.HC))
rob.se.m3.HC <- sqrt(diag(cov.m3.HC))
rob.se.m4.HC <- sqrt(diag(cov.m4.HC))
rob.se.m5.HC <- sqrt(diag(cov.m5.HC))
rob.se.m6.HC <- sqrt(diag(cov.m6.HC))
rob.se.m7.HC <- sqrt(diag(cov.m7.HC))
rob.se.m8.HC <- sqrt(diag(cov.m8.HC))


# We can now run our interaction models from Table E2
m9 <- plm(loan_rate ~ 
            period + log_covid + estab_covid_pct + log_pop + 
            black + hispanic + pov +
            party + hwm + hsb + happ + ssb + pvi_abs +
            log_don + log_estab_large +
            rural_sub + sparse_sub + dense_sub + urban_sub + purban +
            state_AL + state_AR + state_AZ + state_CA + state_CO + 
            state_CT + state_DE + state_FL + state_GA + state_HI + 
            state_IA + state_ID + state_IL + state_IN + state_KS + 
            state_KY + state_LA + state_MA + state_MD + state_ME + 
            state_MI + state_MN + state_MO + state_MS + state_MT + 
            state_NC + state_ND + state_NE + state_NH + state_NJ + 
            state_NM + state_NV + state_NY + state_OH + state_OK + 
            state_OR + state_PA + state_RI + state_SC + state_SD + 
            state_TN + state_TX + state_UT + state_VA + state_VT +
            state_WA + state_WI + state_WV + state_WY +
            period*black,
          data = panel, model = "random")

m10 <- plm(loan_rate ~ 
            period + log_covid + estab_covid_pct + log_pop + 
            black + hispanic + pov +
            party + hwm + hsb + happ + ssb + pvi_abs +
            log_don + log_estab_large +
            rural_sub + sparse_sub + dense_sub + urban_sub + purban +
            state_AL + state_AR + state_AZ + state_CA + state_CO + 
            state_CT + state_DE + state_FL + state_GA + state_HI + 
            state_IA + state_ID + state_IL + state_IN + state_KS + 
            state_KY + state_LA + state_MA + state_MD + state_ME + 
            state_MI + state_MN + state_MO + state_MS + state_MT + 
            state_NC + state_ND + state_NE + state_NH + state_NJ + 
            state_NM + state_NV + state_NY + state_OH + state_OK + 
            state_OR + state_PA + state_RI + state_SC + state_SD + 
            state_TN + state_TX + state_UT + state_VA + state_VT +
            state_WA + state_WI + state_WV + state_WY +
            hispanic*period,
          data = panel, model = "random")

m11 <- plm(loan_rate ~ 
             period + log_covid + estab_covid_pct + log_pop + 
             black + hispanic + pov +
             party + hwm + hsb + happ + ssb + pvi_abs +
             log_don + log_estab_large +
             rural_sub + sparse_sub + dense_sub + urban_sub + purban +
             state_AL + state_AR + state_AZ + state_CA + state_CO + 
             state_CT + state_DE + state_FL + state_GA + state_HI + 
             state_IA + state_ID + state_IL + state_IN + state_KS + 
             state_KY + state_LA + state_MA + state_MD + state_ME + 
             state_MI + state_MN + state_MO + state_MS + state_MT + 
             state_NC + state_ND + state_NE + state_NH + state_NJ + 
             state_NM + state_NV + state_NY + state_OH + state_OK + 
             state_OR + state_PA + state_RI + state_SC + state_SD + 
             state_TN + state_TX + state_UT + state_VA + state_VT +
             state_WA + state_WI + state_WV + state_WY +
             pov*period,
           data = panel, model = "random")

m12 <- plm(loan_rate ~ 
             period + log_covid + estab_covid_pct + log_pop + 
             black + hispanic + pov +
             party + hwm + hsb + happ + ssb + pvi_abs +
             log_don + log_estab_large +
             rural_sub + sparse_sub + dense_sub + urban_sub + purban +
             state_AL + state_AR + state_AZ + state_CA + state_CO + 
             state_CT + state_DE + state_FL + state_GA + state_HI + 
             state_IA + state_ID + state_IL + state_IN + state_KS + 
             state_KY + state_LA + state_MA + state_MD + state_ME + 
             state_MI + state_MN + state_MO + state_MS + state_MT + 
             state_NC + state_ND + state_NE + state_NH + state_NJ + 
             state_NM + state_NV + state_NY + state_OH + state_OK + 
             state_OR + state_PA + state_RI + state_SC + state_SD + 
             state_TN + state_TX + state_UT + state_VA + state_VT +
             state_WA + state_WI + state_WV + state_WY +
             log_estab_large*period,
           data = panel, model = "random")

cov.m9.HC <- vcovHC(m9, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m10.HC <- vcovHC(m10, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m11.HC <- vcovHC(m11, method = "white1", type = "HC3", 
                    cluster = "group")
cov.m12.HC <- vcovHC(m12, method = "white1", type = "HC3", 
                    cluster = "group")

rob.se.m9.HC <- sqrt(diag(cov.m9.HC))
rob.se.m10.HC <- sqrt(diag(cov.m10.HC))
rob.se.m11.HC <- sqrt(diag(cov.m11.HC))
rob.se.m12.HC <- sqrt(diag(cov.m12.HC))

stargazer(m9, m10, m11, m12,
          se=list(rob.se.m9.HC, rob.se.m10.HC,
                  rob.se.m11.HC, rob.se.m12.HC),
          type = "html", out = "PoP.html")


# Now we can generate the different plots.
# First, grid of coefficient plots based on full model m8 (FIGURE 5).

# Color plots for presentation at JFK, Berlin
set_theme(theme_classic())

# Control variables
con.coef.f5 <- plot_model(m8,
      terms = c("log_covid", "estab_covid_pct", "log_pop"),
      title = "",
      vcov.fun = "HC",
      vcov.args = list(type="HC3"),
      show.p = T,
      show.values = T,
      value.offset = 0.2,
      value.size = 4,
      dot.size = 5,
      line.size = 1,
      digits = 3,
      vline.color = "grey")

con.coef.f5 <- con.coef.f5 + 
  theme_classic() + ylim(-.2, .2) +
  scale_x_discrete(labels=list(log_covid = "COVID-19 cases (ln)",
                               estab_covid_pct = "Vulnerable Businesses (ln)",
                               log_pop = "Population (ln)")) +
  ggtitle("Controls") +
  theme(plot.title = element_text(color = "black")) +
  theme(axis.text = element_text(size = 12))

# Political variables
pol.coef.f5 <- plot_model(m8,
      terms = c("party","hwm","hsb","happ","ssb", "pvi_abs"),
      title = "",
      vcov.fun = "HC", 
      vcov.args = list(type="HC3"),
      show.p = T,
      show.values = T,
      value.offset = 0.4,
      value.size = 4,
      dot.size = 5,
      line.size = 1,
      digits = 3,
      vline.color = "grey")

pol.coef.f5 <- pol.coef.f5 + 
  theme_classic() + ylim(-.1, .1) +
  scale_x_discrete(labels=list(party = "Republican",
                               hwm = "House Ways & Means",
                               hsb = "House Small Business",
                               happ = "House Appropriations",
                               ssb = "Senate Small Business",
                               pvi_abs = "PVI (absolute)")) +
  ggtitle("Conventional Politics") +
  theme(plot.title = element_text(color = "black")) +
  theme(axis.text = element_text(size = 12))


# Business power variables
bp.coef.f5 <- plot_model(m8, 
     terms = c("log_don", "log_estab_large"),
     title = "",
     vcov.fun = "HC", 
     vcov.args = list(type="HC3"),
     show.p = T,
     show.values = T,
     value.offset = 0.15,
     value.size = 4,
     dot.size = 5,
     line.size = 1,
     digits = 3,
     vline.color = "grey")

bp.coef.f5 <- bp.coef.f5 + 
  theme_classic() + ylim(-.1, .1) +
  scale_x_discrete(labels=list(log_don = "Campaign Donations (ln)",
                               log_estab_large = "Large Establishments (ln)")) +
  ggtitle("Business Power") +
  theme(plot.title = element_text(color = "black")) +
  theme(axis.text = element_text(size = 12))

# Socioeconomic variables
soc.coef.f5 <- plot_model(m8, 
      terms = c("period", "black", "hispanic", "pov"),
      title = "",
      vcov.fun = "HC", 
      vcov.args = list(type="HC3"),
      show.p = T,
      show.values = T,
      value.offset = 0.2,
      value.size = 4,
      dot.size = 5,
      line.size = 1,
      digits = 3,
      vline.color = "grey")

soc.coef.f5 <- soc.coef.f5 + 
  theme_classic() + ylim(-.5, .5) +
  scale_x_discrete(labels=list(
    period = "Period",
    black = "Black Population",
    hispanic = "Hispanic Population",
    pov = "Poverty Rate"
  )) +
  ggtitle("Socio-economic Context and Time") +
  theme(plot.title = element_text(color = "black")) +
  theme(axis.text = element_text(size = 12))

# Generate grid
grid.arrange(con.coef.f5, pol.coef.f5, bp.coef.f5, soc.coef.f5, nrow=2)

# We used the same commands from sjPlot and ggplot2 to plot the interaction
# models in FIGURE 6; using models m9-m12.

### END OF FILE ###
###########################################################################
###########################################################################
###########################################################################